Optimal stirring strategies 
for passive scalar mixing 



ZHI LIN 1 , JEAN-LUC THIFFEAULT 12 
and CHARLES R. DOERING 13 

institute for Mathematics and Its Applications, University of Minnesota, Minneapolis, MN 

55455 

2 Department of Mathematics, University of Wisconsin, Madison, WI 53706 

3 Department of Mathematics, Department of Physics, and Center for the Study of Complex 
Systems, University of Michigan, Ann Arbor, MI 48109 

We address the challenge of optimal incompressible stirring to mix an initially inhomoge- 
neous distribution of passive tracers. As a quantitative measure of mixing we adopt the 
H^ 1 norm of the scalar fluctuation field, equivalent to the (square-root of the) variance 
of a low-pass filtered image of the tracer concentration field. First we establish that this 
is a useful gauge even in the absence of molecular diffusion: its vanishing as t — > oo is 
evidence of the stirring flow's mixing properties in the sense of ergodic theory. Then we 
derive absolute limits on the total amount of mixing, as a function of time, on a periodic 
spatial domain with a prescribed instantaneous stirring energy or stirring power budget. 
We subsequently determine the flow field that instantaneously maximizes the decay of 
this mixing measure — when such a flow exists. When no such 'steepest descent' flow 
exists (a possible but non-generic situation) we determine the flow that maximizes the 
growth rate of the H~ l norm's decay rate. This local-in-time optimal stirring strategy 
is implemented numerically on a benchmark problem and compared to an optimal con- 
trol approach using a restricted set of flows. Some significant challenges for analysis are 
outlined. 



1. Introduction. 

The enhancement of mixing by stirring in incompressible flows is an important phe- 
nomenon in the natural sciences and in engineering applications. A natural question to 
pose is, how efficient a mixer can an incompressible flow be? This fundamental question, 
more precisely posed, is the subject of this paper. 

In principle, given an appropriate quantitative measure of mixing along with suitable 
constraints on the accessible class of flow fields, the most efficient mixing strategy may 
be determined by solving an optimal control problem. In practice this may be difficult, 
so it is useful to consider other approaches that might more easily be implemented, at 
least theoretically or computationally. Moreover, it is always useful to know absolute 
limits on how fast mixing could ever be achieved subject to the relevant constraints. 
Such bounds provide a scale upon which particular strategies may be evaluated to gauge 
their effectiveness. Here we propose and analyze a theoretical scenario with a particular 
mixing measure subject to particular constraints where these issues can be investigated 
analytically and via direct numerical simulation. As will be seen, there is an interesting 
gap between the best available a priori analysis and the simulation results. 

The rest of this paper is organized as follows. In the next section we introduce the 
basic notions and define the specific problems to be studied. Rigorous bounds on the 
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rate at which mixing measures can decay for constrained stirring flows are derived in 
section 3, and an optimal mixing strategy is described in section 4. In the final section 
5 we report computational implementations of the optimal strategy and discuss open 
challenges suggested by the results. 



2. Problem description. 

Consider the advection of a passive scalar field 9(x, t) by a smooth incompressible flow 
field u(x,t) described by the partial differential equation 

d t 6 + u-V6 = (2.1) 

along with initial condition 9(x, 0) = 0q(x). The stirring field u and the initial distribu- 
tion 8o are periodic in the d-dimensional domain [0, L] d so the subsequent solution 9 is 
as well. Without loss of generality 9q and 9 are spatially mean-zero: 



(9{-,t)) = ±j f 9{x,t)dx = 0. (2.2) 

L ■l[0,L] d 



We also restrict attention to spatially mean-zero flows, i.e., (ui(-, t)) — for i = 1, . . . , d. 

The goal of effective stirring is to redistribute the passive tracer density to achieve 
a maximal rate of mixing as quantified by the relevant mixing measure. To gauge the 
effectiveness of the stirring we adopt as a mixing measure the H~ l norm on mean-zero 
functions where, for any real parameter a, the H~ a norm is a weighted sum of Fourier 
coefficients of the scalar field: 

\\8(;t)\\ 2 H-« = ll|V|-°0(-,t)||i a = $>- 2o |0 fe (i)| 2 (2.3) 

where k = \k\ and 

fc (t) = -L f e- ik *e(x,t)dx. (2.4) 

L 1 J[0,L] d 

The operator |V| _a generally acts in Fourier space as multiplication by k~ a although 
when a is an odd integer V _a is naturally defined as multiplication by —ik/k^- a+1 \ 

We focus in particular on the H^ 1 norm which is related to the large-scale mixing 
measure previously studied by a subset of the authors ( Doering fc Thiffeault|2006 Shaw 



et al. 2007). It measures the variance of a low-pass-filtered image of the concentration 
field: the smaller the H^ 1 norm is, the more homogeneous the scalar field is on large 
spatial scales. In many applications molecular diffusion, implemented mathematically 



by an additional nA9 term on the right hand side of the equation (2.1), dissipates the 
variance, i.e., the H° norm, of the scalar field. But even in the absence of molecular 
diffusion when the H norm is conserved, or when the diffusion is ineffective on the 
length and time scales of interest, the relatively large-scale structures in the scalar field 
may nevertheless decay when a 'mix no rm' like the H^ 1 norm is employed. This idea 



was introduced by Mathew et al. (20051 for the H~ x / 2 norm and is extended to other 
norms here. 

For a > the norms H~ a provide a quantitative measure of mixing in the sense of 
ergodic theory. To see this consider Lagrangian particle trajectories X(t) defined by 

d ^ = u(X(t),t) (2.5) 
with random initial condition X(0) distributed according to the density po(x). Then the 
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tracer particle positions are distributed according to the solution p(x,t) of 

d t p + u ■ Vp = (2.6) 

with initial condition po{x). Incompressible advection conserves the variance of p(x,t): 
the mean-zero field 8 = p — L~ d also satisfies the advection equation, and multiplying 
(2.1) by 9, integrating over the domain, and integrating by parts yields d||0||£2/dt = 0. 



Hence the variance of the density does not measure mixing. Rather, the stirring field 
u(x,t) is called mixing if for every periodic square-integrable function g(x) on [0,L] d , 



lim 

t— ¥00 



[0,L]' 



g(x)p(x,t)dx = (g). 



(2.7) 



See, for example, Lasota & Mackey (1994). The utility of the H a norms (2.3) arc 



indicated by the following theorem, an extension of Mathew et al. ( 2005 ) . 

Theorem: Suppose the spatially mean-zero function 9(x,t) is bounded uniformly in 
L 2 ([0,L} d ) for all t > 0. Then 



lim 



[0,1]' 



g(x)6(x,t)dx = V.g e L 2 



lim 



,t)\\H-* = 0, for any a > 0. 



See the appendix for an elementary proof that applies as well to many other measures 
which could serve just as effectively as a mix norm in this regard. 

The upshot is that more rapid self-averaging characteristic of the intuitive notion of 
mixing is indicated by more rapid decay of the H~ a norm. Alternative measures have also 
been used to characterize mixing and the associated control problem. See, for example, 



Sharma & Gupte 


1997); D'Alessandro et al. ( 1999|); |Vikhansky ( 


2002); 


Schumacher & 


Sreenivasan ( 2003 ) 


; Balogh et al. (20051; Mathew et al. (20051; Thiffeault & Finn (2006); 


Mathew et al. 


(2007); Constantin et al. 


(2008); Cortelezzi et al. 


(2008 


); Thiffeault & 


Pavliotis (2008 


; Gubanov & Cortelezzi ( 


2010). The H 1 norm adopted here, however, 



allows for the development of a particularly straightforward and operational stirring 
strategy. 

Constraints must be imposed upon on the available flow fields in order to formulate an 
optimization problem. We focus on constraints of bounded instantaneous kinetic energy, 
proportional to the L 2 norm of the velocity ||«(-,£)||^ 2 , or bounded instantaneous power 
dissipation in the flow, which for Newtonian fluids with periodic boundary conditions is 
proportional to the H norm of the velocity || Vu(-, t)\\% = ££j=i ||ft«iO,*)Hia- That 
is, we consider flow fields u satisfying either 



[Q,L] d 



|u(a:,i)| 2 da: = U 2 L d 



or 



d 



f |Vtt(x,t)| 2 dx = / i 2^(d iUj 



djUi) 2 dx 



L" 



(2.8) 



(2.9) 



*[0,L] d J[Q,L] a 

defining, respectively, the root-mean-square velocity U — (lii] 2 } 1 / 2 or rate of strain r _1 = 
(iVitl 2 ) 1 / 2 of the stirring. Given either constraint two natural questions are: 

(I) What flow minimizes the mixing measure evaluated at a final time t& a > 0? 

(II) What flow maximizes the instantaneous decay rate of the mixing measure? 

These questions have different answers implying different 'optimal' stirring strategies. 



Mathew et al. (2007) studied problem (I) and solved it numerically for a limited set of 
flow field configurations using the H~ r l 2 norm as the mixing measure. They computed 
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the controls in the form of a time varying linear combination of two simple cellular 
flows. That approach is global in time since it requires keeping track of the complete 
evolutionary history of the system within the interval [0, ffij. 

In this paper we address question (II) and consider flows that produce the steepest 
descent of the H^ 1 mixing measure at each instant in time. This local-in-time strategy 
identifies an optimal mixing flow u(x, t) at time t using only a snapshot of the scalar field 
9(x,t) at that instant. It should be noted that for the optimization problem (I) seeking 
to minimize the mixing measure at a final time, natural constraints might also be the 
total action (the integral of \\u(-, from t = to ta n ) or the total energy (proportional 
to the time integral of || Vit(-, tjW?, from t = to tn n ). 



3. Absolute bounds on mixing rates. 

As a starting point it is useful to identify absolute limits on the rate at which scalar 
fields might be mixed by fluids satisfying the instantaneous energy or power constraints. 
Toward this end we note that incompressible advection conserves not only the variance of 
9 but also that the (weak) maximum principle ensures that the L°° norm (the supremum 
of \9\ over the spatial domain) is conserved as well, i.e., i)||t«> = H^olU 00 a t every 
time t > 0. 

-A~ 1 9, integrat- 



First consider the fixed energy constraint (2.8). Multiplying (2.1 1 by 



2 

\h- 



= -2 



i u 



V(A" 



ing over the spatial domain, and integrating by parts implies 

The Holder and Cauchy-Schwarz inequalities then give 

^||f||ir-i>-2HU 3 || 
and dividing both sides by 2 an d integrating in time yields 

\\0(-,t)\\ H -i > ||0o|k-i U L d ' 2 1 1 9 1 1 l~ t. 



dx. 



-2 UL d / 2 \\9 \\ L ~ \\0\\ H -i 



(3.1) 



(3.2) 



(3.3) 



This rigorous estimate does not rule out perfect mixing as measured by the H 1 norm 
after a finite time, but it does bound the absolute minimum mixing time from below by 



1 



70UH- 1 



ULW \\9 \\ 



2ttU' 



singling out a length scale 



l = 2tt- 



|V-^ | 2 ) 1/2 



?0||L° 



(< L due to the Poincare and Holder inequalities) 



(3.4) 



(3.5) 



characterizing the spatial extent or 'size' of initial inhomogeneities in the scalar field. The 
lower bound in (3.3 1 simply states that under the constant energy constraint, the time it 
takes to achieve complete mixing is at least as long as the time it takes to transport scalar 
inhomogeneities across the characteristic distance ~ Iq at the rms speed U of the flow. 
Whether or not this limiting mixing rate can actually be achieved, or even approached, 
by any suitably constrained stirring flow remains to be seen. 

The analysis is rather different for flows subject to the fixed power constraint (2.9). 
For convenience we define the filtered scalar field 



<p(x,t)= (A-'9) (x,t) 



^2 &ih " 

fe#0 



'k~ 2 9 k (t). 



(3.6) 
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Then starting from (3.1) more integrations by parts yield 

d 



1| 
df 



2 

Iff" 1 



df 



|Vv»| 



2 
L- 



d 



dxidxj 



and the Holder and Cauchy-Schwarz inequalities imply 
-2|M| L oo ||V«|| L 2 IIA^H^ = 



1| 
df 



Iff- 1 



> 



2|MU° 



(3.7) 



(3.8) 



In order to close the differential inequality it is necessary to bound the sup norm 

in terms of the H +1 norm of tp, i.e., || V<£>||/,2 = || V~ 1 9\\ L 2 — \\6Wh- 1 j ancl some conserved 

(or otherwise a priori bounded) quantities. This is possible in 2 and 3 spatial dimensions. 

First consider d — 3. For mean-zero functions on the 3-torus, there exists an 0(1) pure 
number C3 so that 



C 3 



9o||^ 2 . 



For an elementary proof of the inequality see Doering & Gibbon (1995). Thus 



-I 
df 1 



Iff- 1 



> 



C3L 3 / 2 



|0o 



1 3/2 

L 2 



-1/2 



ff- 



and 



> 



^ollff- 1 



1 - 



3C 3 
2r 



L 1 1 6*o 1 1 
ll^ollff- 1 



3/2 



X t 



(3.9) 



(3.10) 



(3.11) 



as long as the term in brackets is non-negative. Again the rigorous analysis does not rule 
out perfect mixing in a finite time but bounds the minimal mixing time from below by 



where 



£ = 2tt^ 



AlUff- 1 

l^olk 2 



2 ( £0 \ 
3C 3 \2itL) 



[< L by Poincare's inequality) 



(3.12) 



(3.13) 



is another length scale characterizing the size of inhomogeneities in the initial distribu- 
tion. 

The linear dependence of t m i X on t is not surprising. This rigorous estimate is, however, 
more than naive dimensional analysis because the lower bound on the mixing time allows 



for a nontrivial i-dependence. The minimal mixing time estimate in (3.12) suggests that 



the availability of larger domain length scales L for the flow to access may in fact facilitate 
mixing, perhaps by allowing for more effective 'folding' to accompany judiciously localized 
'stretching'. The suggestion (which remains a conjecture at this point) is that if N 3 copies 
of the same initial distribution with basic scale £q are assembled into a larger periodic 
domain [0, N x L} 3 , the scalar might possibly be mixed faster by a flow with the same 
rms rate of strain than on the elementary cell [0,X] 3 . 

A similar sort of system size dependence is suggested in two dimensions. The d = 2 



analog of (3.9) is (again, see Doering & Gibbon (1995) for an elementary derivation) 



MIl- < C 2 ||V^|| iM /l + log 



L\\Atp\\ L2 
27r||V(p|| L2 



c 2 \\e\ 



H- 




?0||L 2 



2tt 



Iff- 
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for an 0(1) constant C%. Thus 

C 2 L „„ „ 



-I 



\H~ 1 



> - 




L\\8 \\ L 2 
2ir\W\ 



H- 1 



(3.14) 



This differential inequality also does not prevent from vanishing in finite time, 

but it guarantees that cannot vanish before the absolute minimum mixing time 



tn 



2^02 



log - 



where 1$ is defined in (3.13). For L 3> io the minimum mixing time is 

-1/2 



in 



4 

T X — X 

L 



log 



(3.15) 



(3.16) 



again allowing for the possibility that it may take less time to mix within a larger volume. 

In contrast to these potentially finite-time mixing scenarios, if the flow field is con- 
strained to have a uniformly bounded (in space and time) rate of strain then it can decay 
at most exponentially. In any spatial dimension if || (Vit) sym ||_Loc < 7 < 00 then 

U ' ^^O- (Vm) ■ V _1 0das > -2 7 ||6'||^- 1 (3.17) 



df 



2 

Iff- 1 



and Gronwall's inequality implies an exponential lower bound on the mixing measure: 



||0(-,t)|| H -i > \\64h-^ 



-it 



(3.18) 



Thus if the rate of strain is bounded then the mixing rate cannot increased by increasing 
the domain scale (holding all other constraints fixed). Therefore if there is any real 
precision regarding the L-dependence in the estimates (3.12) or (3.16) then it depends 
crucially on a fixed-power flow's freedom to locally intensify the strain. 



4. Optimal stirring. 



Recalling (3.1) we write 
d 

dt' 



|V-^ = - 2 /»-(^)d a , = -2/„.p(,V^ (4.1) 

where P(-) is the projector onto divergence-free fields defined by 

P(v) =v - VA-\V -v). (4.2) 



Then with either the fixed energy (2.8) or fixed enstrophy (2.9) constraints it is easy to 



see that the optimal mixer maximizing the rate of decay of the mix-norm H 1 is 

"*= V(X)ry/ 2 (4 ' 3) 

or 

_ 1 -A- 1 P(6»V(^) 



~ p r (iV-iP^V^)! 2 ) 1 ^ (4 ' 4) 
if the norm in the denominator does not vanish. So long as the relevant optimal stirring 
exists, the ideal instantaneous stirring strategy is to implement it at each instant of time. 
But if either of the norms in the denominators vanishes then P(0Vy>) = throughout 
the domain and no incompressible flow can instantaneously decrease the mix-norm H —1 . 
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A sufficient (but not to our knowledge necessary) condition for such degeneracies is 
that the scalar field satisfies 

A0 = F(6), (4.5) 

which includes cases where is an eigenfunction of the Laplacian. Such situations arise 
naturally as convenient initial conditions but we suspect that they are non-generic dis- 
tributions among solutions of the advection equation. Nevertheless when and if this 
situation develops some other strategy must be adopted to stir the fluid. 

The natural thing to do when P(0V<^) = is find the flow that maximizes the rate of 
increase of the rate of decrease of the norm. That is, seek the flow that minimizes 

i2 r 

^2^11^-1=2 / [«• V<pV0-tt- (u- V0)A -1 (tt- V0)] dsc. (4.6) 

The optimal incompressible u here solves the eigenvalue problem 

Am =f((u- V8)V(p + (u ■ Vtp)V8 - 2[A~ 1 (« ■ V0)]Vfl) (4.7) 

for fixed energy constraint or 

Am = -A _1 p((u • V8)Vtp + (u ■ Vip)V8 - 2[A~ 1 (w • V0)]V0~) (4.8) 

for fixed power constraint. In either case we seek the eigenvector (field) corresponding to 
the minimum eigenvalue A_ < to use as the stirring field momentarily until P(0V</j) ^ 
0. 

These are generally difficult eigenvalue problems. To make some analytical progress 
consider the special case 8q(x) = sin(fca;) that we will use as an initial condition in the 
computational test reported in the next section in spatial dimension d = 2. Introducing 



the stream function so that u — dip/dy and v = —dip/dx and taking the curl of (4.7) 



- AA?/> = -2V xFv = -2V x (v - VA _1 (V • v)) = -2V x v = 2^ (4.9) 

oy 

where v(x,t) = (it • V0)V<p + (u ■ Vip)V8 - 2[A~ 1 (m • V0)]V0. Writing ip(x,y) = 
n ipmne 27T ^ mx+ny ^ L produces a matrix equation for the Fourier coefficients i/w. 



With the fixed energy constraint (2.8 1 finite-resolution (n, m < N) numerical solutions 
to this matrix equation show that the minimum eigenvalue is always associated with the 
A^ th -harmonic, suggesting that the optimal flow depends on the imposed Fourier cut- 



off. Focusing instead on the fixed enstrophy constraint (2.9|, it is readily shown that to 



leading order in the initial data wavenumber k the minimum eigenvalue we seek is 

mhiA mn (fc) = \ m *,i{k) ~ ~ [1 + O ((kL)- 4 )] (4.10) 

and the associated eigenfunction, modulo an arbitrary phase shift in y, is approximately 
proportional to the first harmonic, i.e., u oc cos(2ny / L) and v = 0. 

Finally, we remark that unless the scalar field is a superposition of the lowest available 
modes it can be unmixed. Indeed, when P(0V</?) ^ then simply reversing the optimal 
flow momentarily unmixes the scalar as defined by the H^ 1 norm, and the same is true 
in the exceptional situations with P(0V<^) = 0. 



5. Computational tests and discussion. 

The first test of the optimal stirring strategy uses initial scalar distribution 8q(x) = 
sinx (with £q/L = 1) in a domain of size L = 2% in d = 2 spatial dimensions. We 
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Figure 1. (a) Evolution of mixing measures for 6o(x) = sin(a;). All norms are rescaled by their 
initial values, (b) Decay of H~ l for different values of £o/L (solid lines). The £o/L — 1 case uses 
8o — sinx, while the Iq/L = 4/5 initial distribution is a sum of 5 low-wavenumber harmonics, 
6> = sin x+ 1.92(-0.542 cos 2x+0.8267 sin y-0.6592 cos 2y+0.3998 sin x cos J/+0.6516 cos x sin y) . 
The £o/L = 2/3 data used 9o(x) = sinx + \/5/3 sin 3y. The £o/L = 1/2 case uses discontinuous 
(non-periodic) initial data, 6o = sin[2(l + 0.5 sin w):r], and the lo/L — 1/3 initial data is sin 3a;. 
The dot-dash lines illustrates the conservation of I? norms with different data in the simulations. 



implemented the fixed power constraint with r 1 = 6.25 x (27r) 2 , which is equivalent 
to the amplitude of the bi-component control utilized by Mathew et al. (2007) . Figure 



1(a) shows the evolution of the scalar mix norms. With the optimal mixing protocol 
proposed here, after an initial transient where the flow is chosen to maximize , 
a robust exponential decay of the H^ 1 norms emerges. In order to compare with the 
two-component flow optimal control results of Mathew et al. { 2007 ) we also computed 
the if -1 / 2 norm that they utilized and reported. Not unexpectedly, expanding the set of 
available flows from two possibilities to everything within the constant power constraint 
allows for faster mixing. Nevertheless the local-in-time optimization strategy employed 
here leaves room for further improvement. The conserved L 2 norm was also monitored 
as a numerical check. 

To check for a possible dependence of the mixing rate on £q/L we also considered 
different initial data, and the results are displayed in Figure 1(b) As is evident, the 



mixing scheme generates robust exponential decay of the mix-norm. Exponential decay 
at a rate independent of the initial data is the sort of dynamics that might be expected 
for a uniform boundedness constraint on the rate of strain, but not the finite-time mixing 
that the best available analysis allows for when only the root-mean-square rate of strain 
is bounded. It is possible that the rigorous analysis may only be saturated by a global-in- 
time, i.e., optimal control, stirring strategy but this discrepancy between rigorous analysis 
and simulation constitutes a major open question. Figure [2] shows the snapshots of the 
scalar field evolution with initial distribution 0q(x) — sinx under the local fixed-power 
optimal mixer. The optimal flow generates a suggestively self-similar cascade of the scalar 
fluctuations to small scales. Other initial conditions generate similar cascades. 
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Appendix: Proof of theorem 

We use the norm ( 2.3 1 for mean-zero functions, \\9\\ H -. = (£ fc kr 2a \9 k \ 2 ) ' \ although 
the proof easily generalizes to other norms. Suppose that 9(-,t) is uniformly bounded 
in L 2 ([0,L] d ), so that \\0(-,t)\\ L2 < C, and lim^^ ||0(-,t)|| H _ a for some a > 0. 
Then for any g E 1? ( [0, L] d ) , 



[0,L] d 



I g dx 



fJk 



< 



k<K 



k>K 



E* 2a 



9k\ 



k<K 



1/2 



L 2 



1 /2 

Given e > 0, first choose K(e) such that (^2k>K(t)\9k\ 2 ) < e /2C, then choose T(e) 
such that ||0(-,T(e))|| ff _ o < \e (J2k<KU) k^fa] 2 )' 1 ' 2 , for t > T(e). Then 



E 

fc>if 



1/2 



[0,L]< 
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< Ui + c- 1 
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e < e, 



t > T(e), 



which implies that 9 converges weakly to zero as t — > oo. (This is true even for a — 0.) 
Conversely, suppose ||0(-,t)||£2 < C for all t and lim^oo J, Q L , d 9gdx — > for all g E 

L 2 ([0, L] d ). By choosing g = cxp(-ik-x) we see that all the Fourier coefficients 6k(t) — > 

as t — y oo. Also, because ||0(-, t)\\\ 2 = Efel#fe(*)| 2 < C 2 > each l^fe(*)l < c for a11 1 - 
Thus 



2 



E fc" 2o i^i 2 + E fc ~ 2a i^i 2 < E fc ~ 2a i^i 2 + ^ 2a ii^ii 
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(A.l) 
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For any e > choose *T(e) > (2C/e) 1 / 2a so that fc~ 2a ||0|| i2 < fc- 2a C < e/2 for fc > if(e). 
(This requires a > 0). Then because for any finite K, J2k<K k~ 2a \9k(t)\ 2 — > as t — >• oo, 
there exists T(e) < oo such that E fc <A'(e) k~ 2 °|0fe(*)| 2 < e/2 for all t > T(e). From ( [ATT I 
we conclude that ||0|| ff - a < e for all t > T(e), which proves the result. 



